function [y0,yx,yxx,yxxx]= projectionStepExpectedControlClosedForm(var0,varx,varxx,varxxx,horizon,model,setup)

% Save output for the core of the model to compute bond prices
thetaTrans = model.theta';
setup.out  = saveStatesControls(thetaTrans(:),setup);

% Storing var0, varx, etc in thetaVar in the same format as theta, i.e. 
% only distinct terms are saved and scaled by scaleX_in_g
nx = setup.nx;
if setup.appOrder == 1  
    thetaVar = zeros(1,1+nx);
elseif setup.appOrder == 2
    thetaVar = zeros(1,1+nx+nx*(nx+1)/2);
elseif setup.appOrder == 3
    thetaVar = zeros(1,1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6);
end
thetaVar(1,1:1+nx) = [var0(1,1)  varx(1,:).*setup.scaleX_in_g'];
if setup.appOrder > 1
    idx = 1+nx;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            thetaVar(1,idx) = varxx(1,i,j)*(setup.scaleX_in_g(i,1)*setup.scaleX_in_g(j,1));
        end
    end
end
if setup.appOrder > 2
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                thetaVar(1,idx) = varxxx(1,alfa1,alfa2,alfa3)*...
                    (setup.scaleX_in_g(alfa1,1)*setup.scaleX_in_g(alfa2,1)*setup.scaleX_in_g(alfa3,1));
            end
        end
    end
end

% Starting values
nx = setup.nx;
if setup.appOrder == 1  
    theta = zeros(horizon,1+nx);
elseif setup.appOrder == 2
    theta = zeros(horizon,1+nx+nx*(nx+1)/2);
elseif setup.appOrder == 3
    theta = zeros(horizon,1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6);
end
setup.dispOn = 0;
for k=1:horizon
    if k == 1
        setup.thetaVar = thetaVar;
    else
        % Updating thetaVar to apply law of iterated expectations,
        % i.e. r(1)_t = E_t[r_t+1] so r(2)_t = E_t[r_t+2] = E_t[r(1)_t+1]
        setup.thetaVar = theta(k-1,:);
    end
    % Settings for the optimizer
    %lambda              = 10;
    %deltaOption         = 1;
    %setup.objectFuncNum = 4;
    %thetaNLS(k,:) = Projection_LMoptimizer(setup.thetaVar(:),setup,...
    %    setup.tolF,setup.MaxIter,lambda,deltaOption,setup.dispOn);
    
    % Input theteVar for ResidualEquationProjExpectedControl is not used
    % here when we only report Evar_cu
    [~,~,Evar_cu] = ResidualEquationProjExpectedControl(setup.thetaVar(:),setup);
    [theta(k,:),resUnitFree,Yfit] = OLSforProj(Evar_cu,setup.xGrid,setup.appOrder,setup.scaleX_in_g);
    %res = resUnitFree.*Yfit' = Y-Yfit';
    %res = resUnitFree.*Yfit'
    %res'*res
    %theta(k,:)-thetaNLS(k,:)
end

% Undo the scaling
scaleX_in_g    = setup.scaleX_in_g;
y0             = theta(:,1);
yx             = theta(:,1+1:1+nx)./repmat(scaleX_in_g',horizon,1);
yxx            = zeros(horizon,nx,nx);
yxxx           = zeros(horizon,nx,nx,nx);
if setup.appOrder > 1
    idx   = 0;
    for i=1:nx
        for j=1:nx
            if j>=i
                idx = idx + 1;
                yxx(:,i,j)    = theta(:,1+nx+idx)/(scaleX_in_g(i,1)*scaleX_in_g(j,1));
                yxx(:,j,i)    = yxx(:,i,j);
            end
        end
    end
end
if setup.appOrder > 2
    idx   = 0;
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                yxxx(:,alfa1,alfa2,alfa3) = theta(:,1+nx+nx*(nx+1)/2+idx)...
                    /(scaleX_in_g(alfa1,1)*scaleX_in_g(alfa2,1)*scaleX_in_g(alfa3,1));
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %yxxx(:,alfa1,alfa1,alfa3)= yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa1,alfa3,alfa1) = yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa3,alfa1,alfa1) = yxxx(:,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %yxxx(:,alfa1,alfa2,alfa2)= yxxx(:,alfa1,alfa2,alfa3);                
                    yxxx(:,alfa2,alfa1,alfa2) = yxxx(:,alfa1,alfa2,alfa3);                
                    yxxx(:,alfa2,alfa2,alfa1) = yxxx(:,alfa1,alfa2,alfa3);                                
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %yxxx(:,alfa1,alfa2,alfa3) = yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa1,alfa3,alfa2) = yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa3,alfa1,alfa2) = yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa3,alfa2,alfa1) = yxxx(:,alfa1,alfa2,alfa3);
                    yxxx(:,alfa2,alfa3,alfa1) = yxxx(:,alfa1,alfa2,alfa3); 
                    yxxx(:,alfa2,alfa1,alfa3) = yxxx(:,alfa1,alfa2,alfa3);                 
                end                
            end
        end
    end
end



end

